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ABSTRACT 

We analyze a system of parabolic nonlinear equations that describe the 
diffusion of a fully collisional plasma across a strong magnetic field. We 
demonstrate that the solution to this system tends to a time asymptotic 
state which is of space-time separable form, <f>(t)f(x). Furthermore, f(x) 
is independent of the initial conditions and <J> ( t ) depends slightly on the 
initial conditions. The rate of decay of the temporal part is governed by 
a nonlinear eigenvalue problem. Since the equations are considered in a 
bounded domain we are able to analyze the effect of boundary conditions on 
the evolution of the system. Additional effects as radiation, heating, and 
particle injection can also be accounted for. Essential differences 
between the behavior of a fully-coupled system and a scalar equation are 
observed. 
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Derivation of Equations 

We consider the problem of radial diffusion and transport of a 
magnetically confined fully collisional plasma. We are interested in the 
long term behavior of the solution subject to boundary conditions and 
various forcing functions. We verify numerically that the long time 
asymptotics are actually achieved within several transit times. Thus, the 
solution quickly evolves into a universal pattern independent of initial 
conditions . 

In order to simplify the equations, we consider the case that the ions 
and electrons have the same temperature. In this case the mass and energy 
diffusion tensor terms are essentially the one-temperature Braginskii 
equations. Rosenbluth and Kaufman^ present some special solutions for the 
equations in an infinite domain. In this study we shall, in contrast, only 
consider bounded domains. Indeed, the effect of boundary conditions is 
very important. We ignore topological changes of the plasma and assume 

that the plasma extends to the wall. 

The resultant equations are 

n+V*(nu,) = S (la) 

t dp 

P + -g- V-(pu d ) = | V* (KVicT) + Sg. (lb) 

We consider the equations in either a Cartesian or cylindrical frame of 
reference. We then igonore all spatial derivatives except in the x or 
radial direction. The diffusion drift velocity u^ is given by 
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where n is the total (ion plus electron) particle density number, kT is the 
temperature in ev, is the perpendicular resistivity and B is the 

magnetic field which is assumed to be constant and in the z direction. 

We first consider these equations in slab symmetry. We shall see 
lacer that the changes introduced by considering cylindrical symmetry are 
minor. Temporarily ignoring the source terms Sp, Sg and using normalized 
units, (1) can be written as 
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where 



K n 




a 



(2c) 


It follows from (2c) that a and K 1 both depend on the ratio of ion mass to 
electron mass given by Kq. We shall consider two specific cases: 
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Hydrogen: - 60.6, = 20.78, a = 0.888 (3a) 

Deuterium: = 86, » 29.25, a * 0.920. (3b) 

The system (2) is augmented by initial conditions 


n(x,0) = n Q (x), p(x,0) = p Q (x). 


(4) 


We consider (2) in the interval 0 < x < 1. At x = 0 we impose symmetry 
conditions while at the edge x = 1 we impose convective boundary 

conditions. Thus, we have 


n + h. n = 0 
x 1 


p x + h 2 P = ° 


at x = 1, 


(5a) 


p=n=0 at x = 0. (5b) 

x x 


The boundary conditions at the edge (5a) simulate a plasma-wall 
interaction. It follows from (5) that the temperature T at the wall 
satisfies T^ = (h^ - h^T. Hence, hj > h 2 implies heat deposition while 
h^ < h 2 implies heat injection at the plasma edge. 

We can also rewrite (4) as a quasilinear system of equations 


3 
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( 6 ) 



where 
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and 
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-2aK 
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We note that the eigenvalues of A are constant and strictly positive. From 
(3) we see that is a relatively large number between 20 and 30. 
Considering the case of Kj large, the eigenvalues of A are 


a l = 7 K 1 + 0(1) 


and 


a 2 = J + OCl/Kj). 


We conjecture that the long term behavior of the system (2) - (3) is 
given by a time separable form independent of the initial conditions (4). 
Thus, we conjecture that 

lim n(x,t) = ij>^(t)N(x) 

t+00 

(7) 


lira p ( x , t ) = <j) (t)P(x). 
t-*-°° 


Numerical experiments indicate that these asymtotic states are reached 
quite quickly. For single nonlinear diffusion equations, this asymptotic 
behavior is well known. ^ It is also valid for a 2x2 system when the 
diffusion tensor is diagonal. In this study we indicate that it is still 
valid for a 2x2 system with a nondiagonal tensor term. 
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2. Ly panov Functional 

In this section we construct a Lypanov functional for (2) - (5) which 
is a strictly decreasing function of time. Hence, in the absence of 
forcing functions n(x,t) must approach zero. In some cases we can simplify 
the algebra by assuming that is large; however, the results are true for 
all K x > 0. 

Let 


then 


I -I/Kjvo 
Q * In P I > 
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dt J 
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By the Cauchy-Schwarz inequality we have that 


(t)“ 


e 2 1 2 

p < -s- n + -z— p 
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We choose 


e 



1 + a + 


cr( 3+a) 




> 0 . 


It then follows that 
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(9) 


We, therefore, require that 


2a 


( 1 + !) 


a + 7 + (3+a)a/K 


1 


Since by (2c), a = 1 - 7/3K^, we require that 


a > 


8 - 7/3K 


1 


16 - 10/K 1 + 7/3 


^ 1 + _1_ + 0(kJ>. 


( 10 ) 


We next treat the boundary terms using (5). We then conclude that 


W Q 


dx < - 


„ / _ „ A 2 1 2, 

2a J D i Q T n x + — I Px dx 

T> 1 n^ X 6p^ X 


2a( 1 + a)D^(l)Q(l)h 1 


(ID 


where 


Q = 


' n \a 

“T7k7I 5 


a > 0, 


v(f ) 1/2 


A = 2a 1 + x - 


a + 7 + (3+a)a/K^ 


> 0, 


( 12 ) 
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and a also satisfies the inequalities (10). 

Note: 

(1) It is not necessary to choose the exponent of p to be -1/Kj in 
the definition of Q. In fact, one can choose the exponent so that the 
coefficient of p x in (7) is exactly zero. For large this gives an 
exponent of the order of -1/K^ + 0(1/K^). 

(2) The boundary term depends only on h^ and not on i.e., the 
decay of Q is affected by mass deposition but not by pressure boundary 
conditions. If the exponent of p in Q is not exactly -l/K^, then h 2 does 
enter into the energy inequality but with a coefficient of order 1/K^. 

We further note that the right-hand side of (11) depends on spatial 
derivatives of n and p. Hence, as long as n and p are not constant, the 
average of Q must approach zero m the absence of source terms. Since a is 
arbitrary, subject only to the inequality (10), and the exponent of p need 
not be exactly -1/K^ (see Note 1), we conclude that n and p must decay to 
zero. 


3. Decay of Solution 

In the previous section we were able to get a priori estimates on the 
decay of the functional Q, (8), (11). In order to get further information 
on the decay rates, we now assume that the time-separable conjecture (6) is 
true. More exactly we set 
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n(x,t) = c}> 1 (t)N(x) 


p(x,t) = cf> 0 (t)P(x) 


and substitute this hypothesis into (2). As usual we introduce two 
separation constants A^, A 2 « We then find that 


*2.5 d, _1/ 2 

l l 4>1 *2 


X 2 <(> 1* 5 ^ 2 2 


and also 


31 D 1 p ' 1/3 3? ' NP ‘ /3 > + h N * 0 


K 1 3x D 1 nC1 37 (PN t * ) + 7 > 2 p " °> 


with boundary conditions 


^ + h N = 0 
dx 1 


H + h 2 P * 0 


at x = 1, 


( 16a) 


dN _ dP 
dx dx * 


at x = 0, 


(16b) 


where D 1 = N 3 ^ 2 P -1 ^ 2 . 

The solution to (14) is given by 



-9- 


<(>2 = 4>q <t>^> = constant (17a) 

and 

<|>j = (t Q + X* ait) lj/ “, t Q = constant (17b) 

where 

n = X 2 /X 1 , 0) = (3 - n)/2, X* = X^ . (17c) 

It follows that the temporal behavior depends on the ratio of the two 
separation constants. It is easily seen that if one changes scales then 
the values of Xj and X 2 change but the ratio Q is invariant. Hence, Q is 
essentially a nonlinear eigenvalue. The value of 9 is determined by the 
global existence of a solution. In practice 9 is found by numerically 
integrating the system (2) and then calculating the rate of decay. 

It is shown in reference 3 there are four critical values for 

n 

9=3, 7/6, 1,0. 

For S2 > 3 it follows that m < 0 and the diffusion process terminates in a 
finite time. For SI = 3 we have exponential decay while for 9 < 3 we have 
algebraic decay. The value = 7/6 is a bifurcation between processes for 
which heat is added to the system through the boundary. Whenever h 2 > h^ 
then heat enters the system and 9 > 7/6; conversely, when h 2 < hj then heat 
leaves the system and 9 < 7/6. Furthermore, if 1 < 9 < 7/6 then all the 
state variables decay and the temperature has an inverted profile. When 
0 < ft < 1 then there is a thermal instability, and the temperature grows in 
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time even though the pressure and density decay. Finally, when < 0, then 
both the temperature and the pressure grow while the density decays. 


4. Computational Results 

In this section we present some numerical results for the system 
(2). All the computations were performed using an explicit four-stage 
Runge-Kutta type scheme described in reference 3. For the all the graphs 
used we have defined Kq as appropriate for hydrogen (3a). Results for 
deuterium (3b) have been similar. Having fixed Kq the main free parameters 
are h^ and that appear in the boundary conditions at the edge (5a). 

In Figure 1 we present a case in cylindrical coordinates with h^ = 10, 
h 2 = 20, and so hj < h 2 « A sequence of five times are plotted for density, 
pressure, and temperature. The plots are at time steps 1, 100, 250, 500, 
900. The initial profile is linear for the density and exponential for the 
pressure. Other runs indicate that the asymptotic state is independent of 
the initial profile. We see that the solution changes very fast near the 
plasma edge. After 1000 time steps which correspond to a time of about .05 
the asymptotic state is almost reached. In these plots the quantities are 
normalized to be 1 at the center, x = 0. Hence, we only see the spatial 
portion of the solution and have normalized out the decaying temporal 
portion of the solution. We see that the approach of the pressure to the 
steady state is not monotonic. 

In Figure 2 we have increased h^ to 21 without changing 1^ so that now 
hj > h 2 « Comparing Figures 1 and 2 we see that the pressure now approaches 
its asymptotic form more slowly than in the previous case. The final 
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profiles for the density and pressure look similar to Figure 1, but now the 
temperature has an inverted profile. It is interesting to note that the 
temperature approaches its asymptotic form more rapidly than before. We 
have also run this case in Cartesian coordinates and have included 
radiation effects and the graphs look similar to Figure 2. 
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Figure 2a. Normalized plots of density, pressure, and temperature with 
hj = 21, h2 = 20 in axisymmetric cylindrical coordinates. 

Same time steps as in Figure 1. 
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Figure 2b 
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